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ABSTRACT 

A  new  non  parametric  method  for  estimating  the 
locations  and  power  spectral  densities  of  multiple  wide¬ 
band  sources  from  measurements  provided  by  an  array 
of  sensors,  is  described.  The  proposed  method  treats 
with  the  same  case  multiple  sources  and  multipath  pro-  . 
pagation.  Direction-of-arrival  estimation  for  wide-band 
sources  is  obtained  as  a  limiting  case  of  the  proposed 
method. 


L  INTRODUCTION 

In  radar,  sonar  and  seismology  one  is  frequently 
interested  in  determining  the  locations  and  the  spectral 
densities  ("signatures")  of  radiating  sources  from  meas¬ 
urements  provided  by  an  array  of  sensors.  The  signals 
received  by  the  sensors  consist,  in  the  simplest  case,  of 
sources  to  the  sensors,  may  exist. 

The  conventional  method  for  estimating  the  location 
of  the  source  is  based  on  a  two-step  procedure.  First, 
the  time-difference-of-arrival  (TDOA)  of  the  the  propagat¬ 
ing  signal  to  the  different  sensors  are  Estimated  using  a 
generalized  correlator  (see  e.g..  Knapp  and  Carter 
(1976)).  Then,  these  TDOAs  estimates  are  used  to  derive 
the  corresponding  lines-of-position  whose  “intersection" 
yields  the  location  of  the  source.  This  method  has  seri¬ 
ous  shortcomings,  the  major  one  being  the  inability  of 
the  generalized  correlator  estimator  to  cope  effectively 
with  multiple  sources  and  with  multipath  propagation 
(see  e.g..  Carter  ( 1991)). 

An  interesting  attempt  to  overcome  this  problem 
was  outlined  by  V.orf  et  al.  (1979).  Their  basic  idea  was 
to  use  parametric  models  for  the  sources  and  the  addi¬ 
tive  noises  and  then  to  use  system  identification  tech¬ 
niques  to  estimate  these  parameters  from  the  received 
signals.  The  TOOA  estimates  are  then  extracted  from 
these  parameters.  This  idea  was  further  elaborated  and 
developed  by  Porat  and  Friedlander  (1991)  and  Nehorai 
and  Morf  (1992).  The  shortcoming  of  this  method,  ns 
with  any  parametric  method,  is  its  sensitivity  to  the 
assumed  model. 

The  maximum  likelihood  processor  for  estimating 
the  location  of  multiple  wideband  sources '  has  been 
recently  presented  by  Wax  and  Knilath  (1392a).  This 
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processor  requires  the  knowledge  of  the  spectral  density 
matrix  of  the  sources  and  hence  its  applicability  is  lim¬ 
ited  to  the  esses  in  which  such  knowledge  is  available. 

In  the  special  case  that  the  sources  are  narrowband 
and  are  located  in  the  far-held  of  the  array,  the  problem 
degenerates  to  the  estimation  of  the  direction-of-arrival 
and  the  center-frequency  of  the  radisting  sources, 
namely  to  the  2-D  harmonic  retrieval  problem.  This 
problem  has  been  addressed  by  Yfax  et  al.  (1392b).  where 
a  suboptiir.a!  method,  based  on  the  eigenstruclure  of  the 
covariance  matrix  of  the  received  signals,  has  been 
presented.  This  method  is  an  extension  of  the  method 
presented  by  Schmidt  (1979)  (sec  also  Schmidt  (1901)) 
for  the  special  case  that  all  the  sources  arc  co- 
frequency.  A  similar  method  for  this  special  rose,  based 
on  the  eigenstruclure  of  the  spectral  density  matrix,  has 
been  presented  by  Bienvenu  (1973)  (see  also  Bienvenu 
and  Kepp  (1353).  (1991)). 

In  this  peper  a  new  approach  is  presented  'or  the 
problem  of  estimating  the  locations  and  the  spectral 
density  matrix  of  wideband  sources.  The  .'t-pro.-ieh  is 
nor.pe.-arr.etrlc  anc  hence  robust  and  can  cope  with  mul¬ 
tiple  sources  ar.d  multipath  propagation,  i!  is  a  one-slep 
procedure  besed  on  the  estimation  of  the  eigenstruclure 
of  the  spectre!  density  matrix  of  the  received  signals.  !l 
extends  the  approach  of  Schmidt  (1979)  and  Bienvenu 
(1979)  for  the  case  of  wideband  sources  located  in  the 
ncar-ftt'.d  of  the  array.  It  is  shown  that  tiie  eigenvector 
subspece  corresponding  to  the  repeated  sma.'ie -t  eigen¬ 
value  of  the  spsctral  density  matrix  coritshis  all  the 
information  on  the  locations  of  the  impinging  sources. 
Algorithms  for  extracting  this  information  that  enable 
trade-off*  between  resolution  and  accuracy,  ore 
presented.  Simulation  results  that  demonstrate  the  per¬ 
formance  of  the  proposed  algorithms  are  aisc  presented. 


n.  PRC3LKV  FORMULATION 

Assume  that  we  have  m  sensors  and  4  (4  <m) 
sources  distributed  in  the  ptane.  l.et  (x,.y,)  denote  the 
coordinates  of  the  i-th  source.  F.ach  sourc-  >s  assumed 
to  emit  a  zero-mean  widc-sensc-stalionarv  -ignat  that 
propagates  radially  with  speed  e.  Denoting  the  received 
signal  at  the  i-th  sensor  by  r4(f  ).  we  can  write 

r»(0  =  V  aw*»(f-T»)  ♦  n,(f)  (1) 

I*! 

r  t 

-  &  -5-,  l*<  *m 

*  JC 

where 

i«(f  )  a.  the  signal  radiating  by  the  k  -th  source 
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r*  =  the  propagation  time  from  the  Jfc-th 
source  to  the  <-th  sensor 
On  b  the  attenuation  from  the  /t-th  source  to 
the  t-th  sensor 


n,(i)  =  the  additive  noise  at  the  {*tb  tensor 


We  assume  that  the  noises  Jn^f )[«"!,  are  wide- 
sense-stationary  with  zero-mean  and  identical  spectral 
densities,  and  that  they  are  uneorrclated  with  each 

other  and  with  each  of  the  signals  ,  k  =  1 . dj. 

7*  7* 

Since  the  "observation-  interval  T,  =  [-— .  is 

finite,  we  can  represent  the  received  signals  (1)  by  either 
a  Fourier  aeries 


r,(0  =  <2) 


(which  implies  the  periodicity  of  r(  -  )  outside  the  inter¬ 
val  f,)  or  by  the  inverse  Fourier-transform  of  a  Whit¬ 
taker  (sampling)  series: 


r«(t)  =  r-« 


T,/*  £  f?*(u*) 


(«»| — «*) 


(3) 


(which  implies  that  our  processes  vanish  almost  every¬ 
where  outside  the  observation  interval).  In  each  case, 
the  Fourier  coefficients  \Rt{un)  .  s>»  £fl,  where  B  is 
the  bandwidth  of  the  processes  {.  are  given  by 


n<0 

where 

=  — - .  n  tO.  U  •••  (4b) 

Note  that  since  r4(  ■  )  is  real. 

A,(o-*)  =  HiH)  for  all  n  (4c) 

where  •  denotes  the  complex  conjugate.  Thus  we  need 
only  to  consider  positive  “frequencies",  i.e..  u,  with 
n  >  0. 

Taking  the  Fourier  coefficients  of  both  sides  of  (1). 
assuming  tmOCF.  we  obtain 

/ftC'j.)  *  £  ot*5a(e»ll)«'iv«  ♦  *»(«*,)  <  =  (5) 

*•» 

or  rewritten  in  matrix  notation 

R(d„)  *  A(o,)S(u„)  +  N(s>»)  (6a) 

where  R(s»„)  .  S(u„)  and  N(u«)  are  the  mxt  vec¬ 
tors 


*.<**.) 

s,M 

*<«J  e 

S(«a)  = 

• 

*m(un) 

*.<««) 

*<««>» 

M«*>] 

land  A(u«)  is  the  rnxrf  matrix 


(6b) 


C 


A(o»)  = 


o,,e 


aMe 


-IVw 


-/svi 


(6c) 


Note  that  the  i-th  column  of  A(uK)  corresponds  to  the 
attenuations  and  delays  of  the  i-th  source  to  the 
different  sensors.  Motivated  by  this  we  define  the  "foe a- 
Ken  vector"  of  the  i-th  source  as. 


“it* 


-f‘VT«t 


-<Vw 


(6d) 


Observe  that  since  S(w»)  and  N(y*)  are  zero- 
mean,  so  is  R(y„).  Thus  it  follows  from  (6)  that  the 
covariance  matrix  of  R(vn)  is  given  by 

£:R(«„)R*(u»)  =  A{»»)£S(aa)8»(«w)A*(a1.) 

trKfa.JK^)  (7) 

where  +  denotes  the  complex-conjugate  transpose. 
Now,  it  is  well  known  (see.  e.g..  Whalen  (1971).  p.  81)  that 
if  the  observation  interval  is  targe  compared  to  the 
correlation  time  of  the  processes  t*‘t(f)K5|.  l**(f)l*»i 
and  t»i(f}|t!i.  then 

£R(s>»)R‘(iv)  *  K (t>n)  (8.a) 

*S(Ma)S»(&>a)  »  P(u„)  (8.b) 

£N(y„)N‘(y*)  *  Q(t>«)  (8.c) 

where  K(y„).  P(a»)  and  Q(un)  are  the  power  spectral 
density  matrices  of  the  processes  Jr4(  •  )J.  {«*(  •){  and 
Jn «(  •  ){.  respectively,  at  the  frequency  «».  Now.  since 
we  have  assumed  that  the  noises  jiqt  -  )|  are  uncorre¬ 
lated  processes  with  the  «ame  spectral  densitirs  it  fol¬ 
lows  that  the  noise  power  spectral  matrix  is  given  by 

9(yn)  =  e*(t>»)l  (9) 

where  c*(s>„)  Is  the  (scalar)  power  spectral  oensity  of 
each  of  the  noises  at  the  frequency  and  1 

is  the  identity  matrix. 

Equation  (7)  can  then  be  rewritten  as 

K(«„)  =  A(--„)P(«JA*(«J  4  c»(ujl  .  (10) 

which  is  the  basic  relation  in  the  forthcoming  analysis. 


ID.  EIGEN  DECOMPOSITION  OKTHESPECTHAl. 

DENSITY  MATRIX 

Observing  the  structure  of  the  spectre!  density 
matrix,  as  given  by  (10).  assuming  that  the  rank  of 
A(»,|)  is  d  (Le..  that  the  "location  vectors”  of  the  d 
sources  are  linearly  independent  )  and  that  P(y„)  is 
positive  definite,  it  can  easily  be  shown  (sec  e.g.. 
Schmidt  (1979)  or  *henvenu  (1979))  that  the  minimal 
eigenvalue  of  K(y»)  is  equal  to  e*(yK)  with  mulliplicity 
m-d,  and  that  the  corresponding  eigenvector  subspace 
is  orthogonal  to  the  columns  of  the  matrix  A («»). 
namely,  to  the  "taenfisn  vectors"  of  the  sources. 

Thus,  denoting  by  J\(y«)|  and  JV« (i-,)|  the  eigen¬ 
values  and  eigenvectors,  respectively,  of  K(va).  where 

A,(»„)  «  X,(u„)  *  «  X»(y») 

it  follows  that 
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=  ••••  =  Kk-*M  =  (h.b) 

and 

lV,(t»B).  i  =  I . m-di  ±  i  =  I . d I .  (11. b) 

where  J-  denotes  orthogonality. 

_  In  practice  K(un)  is  not  known,  so  that  an  estimatt 
K(«w)  must  be  formed,  front  Hie  dgla,._Ihe.  eigenvalue* 
and  eigenvectors  of  K(db)  obey  only  asymptotically  the 
relations  (11).  The  multiplicity  of  the  smallest  m  —d 
eigenvalues  of  K(uB)  (11a)  is  reflected  in  the  estimate 
K(o»)  as  a  “cluster"  of  the  smallest  m  —  d  eigenvalues, 
and  the  orthogonality  condition  (ll.b)  is  reflected  as 

l^si.tt(On)vrM!t  *  o  (12) 

Since  (12)  holds  for  every  «„  €  B.  two  reasonable 
“measures"  of  the  "closeness  to  orthogonality"  over  the 
whole  bandwidth  B  are  either  the  sum  or  the  product  of 
the  individual  "measures"  (12)  at  each  frequency  bin.  in 
the  sense  that 

E  Y  lA*,.*t(oB)V,*(«B)!1  at  0  (13a) 

i-i 

and  also 

n  lA,..*(«a)V, »(«*)!*  »  0  (13b) 

■1.(1  («l 


IV.  DETECTION  OF  THE  NUMBER  OF  SOURCES 

Determining  the  m  -  d  eigenvectors  that 
correspond  to  the  repeated  smallest  eigenvalue  requires 
the  knowledge  of  the  number  of  sources  d.  In  practice 
d  is  usually  unknown  and  hence  it  must  also  be 
estimated  from  the  data. 

The  problem  or  estimating  d  is  equivalent  to  the 
problem  of  determining  the  multiplicity  of  the  smallest 
eigenvalue  of  the  spectral  density  matrix  K(cj„).  This 
problem  can  be  formalized  as  a  hypothesis  test 

Ht(uK)  :  \,(un)  =  =Xm_-(aB) 

A  5  A|(cjb)  /  5  Xm.f(s>n) 


It  can  be  shown  (sec  e.g  Priestley  et  al.  (1973)  )  that  the 
likelihood-ratio  statistic  for  this  problem  is  given  by  the 
ratio  of  the  geometric  mean  and,  the  arithmetic  mean  of 
m  -  d  smallest  eigenvalues  of  K(uB). 


/*("»)  = 


n  A,^)'"-* 

<■1 _ 


(14a) 


Under  the  null  hypothesis  the  statistic  -  2  In  l«(u.)  is 
asymptotically  distributed  as  g*  with  (m  -  d  )*  -  1 
degrees  of  freedom  (see  Priestley  ct  al.  (1973)). 

To  minimize  the  probability  of  error,  this  test 
should  be  performed  for  every  un  C  B.  Xow,  since  the 
Fouricr-cocft'.cient  corresponding  to  diflcrent  frequen¬ 
cies  are  independent  it  follows  that  likelihood  ratio 
statistic  for  the  different  frequencies  are  independent. 
Thus,  the  composite  likelihood  ratio  statistics  for  deter¬ 
mining  the  number  of  sources  is  given  by 

/„  «  UM  (Mb) 

*jiJ  is  distributed  as  g*  with  U[{m  -  d)*  -  1)  degrees 


of  freedom.  The  way  to  implement  this  test  is  to  apply  it 

sequentially  to  d  =  0.1 . m  —  1  and  to  choose  d  as  the 

value  of  d  for  which  —2  In  crosses  the  y* 
significance  level  threshold. 


V.  SOURCE  LOCATION  ESTIMATION 

1*®*  Ar>(i)„)  denote  the  “location  vector" 

corresponding  to  a  source  at  the  point  (z.y)  and  let 
|A*.»(-!i)iVllij  denote  the  collection  of  the  “location 
vectors"  of  ell  possible  locations  j(x.  y){  in  the  plane. 
These  vectors  can  either  be  measured  in  the  field  or 
computed  analytically,  using  the  appropriate  propaga¬ 
tion  model,  if  measuring  is  not  feasible.  The  estimate  of 
the  sources  locations  is  obtained  by  computing  end  plot¬ 
ting 

"l;1  ~  (*'*•’) 

£  £  !  A*.*  (s>B)V|,(yB)|* 

v.ta  »=i 

or 

_____ - 1 - - -  (15b) 

II  E  lAB.»(oB)Y1>(yB)|* 

B  1*1 

as  a  2-D  function  of  (x.y)  .  The  d  peaks  of  these  2-D 
funcUon  constitute  a  good  estimate  of  the  locations 
I  (*!•  y.)  !»*•:•  since  for  this  values  the  denominators  of 
(15)  should  be  "close"  to  zero  according  to  (13) .  The  two 
estimators  give  in  (15)  difier  in  their  resolution  and  sta¬ 
bility  properties.  For  the  first  estimator  (15a).  H  is  elcur 
there  be  a  peak  in  the  point  (x.y)  if  end  only  if 
the  "location  vector"  A,. *(!>„)  corresponding  to  this 
point  is  "closely  orthogonal"  .  at  all  frequencies  £ H  , 

to  the  eigenvectors  }V,(yB).l  =  1 . m-d|  .  The  second 

estimator  (15b),  will  have  a  peak  at  a  point  (z.  y)  even 
if  the  "Iccetion  vector"  A,. „(;>„)  corresponding  to  this 
point  is  closely  orthogonal".  at  only  a  subset  of  lj»:e  fre¬ 
quencies  -nzB  .  to  the  eigenvectors  }V|(o„). 

4  =  1 . m.—d |  .  Thus  the  first  estimator  will  show 

lower  resolution  but  higher  stability  then  the  second 
estimator.  Other  estimators  that  enable  a  more  deiicnlc 
trade-o*  between  resolution  and  stability  ere  discussed 
by  Wax  tt  al.(  1332b). 


VI.  SOURCE  SPECTRAL  DENSITY  ESTIMATION 

Having  at  hand  the  estimates  {(**.  of  Ihc 

sources' location  we  can  construct  an  estimate  A(vn)  of 
the  matrix  A (cB)  simply  by  stacking  together  ihc  i 
"location  vectors"  corresponding  to  the  estimated  posi¬ 
tions  of  the  sf  sources.  Thus 


AC*n)  =  tA;t£((i:B)  A.-<4<(^)) 


(16) 


Xow,  since  the  eigenvalues  and  eigenvectors  of 
A(an)P:i,)A*(aB)  are  MaJ  -  o*(i.B).  Y,(yJ  . 
*  “  m~ .  »i|  then  from  the  well  known  spectral 
theorem,  of  matrix  theory  (see  e.g.,  Strang  (iV'Jij)  on 
estimate  of  A(iiB)P(yB)A*(i)B)  can  be  constructed  by 

A(yn)P?L,)A-(uB)  = 

-5*(“w)]vt(«B)v,*(«B) 

•hor*  *',)  is  the  estimate  of  o*(-B)  •  given  by 


(IVa) 


5’(a»)  =  -rrrit  KM 


(1Tb) 


( 


m 


* 


Rewriting  this  in  matrix  notation  we  obtain 
A(o»)f’(u»)A  *(<*.) 

*  V(«»KA{«»)  -®*(4'a)lIV*(«»)  0*a) 

where 

V(««)  =  •••  VM(oJ]  (18b) 


and 


k»n) 


A«»  ) 


(15c) 


I  ‘  *.(«.)I 

Thus.  solving  (10)  (or  the  desired  spectral  density  matrix 
P(eO .  w*  obtain 


P(«J  =  tA*(«a)M«a))*‘A*(a a)V(a»)[A(aa)  -  **<01} 

•  VH««>A(«w)^A♦(t,,)A<aw)J-,  (15) 


VIL  DDtECTION-OK-ARRIVAL  ESTIMATION  OF  WIDE-BAND 
SOURCES 

Direction-or-arrival  estimation  is  a  limiting  case  of 
source  location  estimation  corresponding  to  sources  at 
~inflnity".  However,  instead  of  deriving  the  direction-of- 
arrival  estimation  algorithm  by  limiting  arguments  let  us 
go  back  to  basics  and  see  how  to  modify  equation  ( 1)  for 
this  case.  Clearly  |rtt|  and  (a*!  are  meaningless  for 
sources  at  infinity.  A  way  around  this  problem  is  to 
define  the  reference  not  at  sources'  locations  but  at 
some  arbitrary  point  in  the  plane,  say.  the  origin  of  the 
coordinate  axes.  Equation  (l)  can  then  be  rewritten  as 

rt(f)  =  t  «‘t (*.)»*(*  -  m(O)  +  «i<0  (20) 

Isef  lsl<n 
2  2 

where 

T\(d»)  E  the  propagating  delay  between  the  4-th  sensor 
and  the  reference  point  for  a  wavefront 
impinging  from  direction  Ok 
n (j, )  =  amplitude  response  of  the  i-th  sensor  to  a 
wavefront  impinging  from  direction  dk 

The  analysis  parallels  that  of  the  general  case  dis¬ 
cussed  in  the  previous  sections,  with  the  only  difference 
that  the  location  vector  degenerates  to  a  direction  vec¬ 
tor  given  by 


Ae,(s>»)  * 


The  direction-of-srrival  estimation  algorithm  proceeds  in 
the  same  manner  as  in  the  general  source  location  prob¬ 
lem.  except  that  (13)  are  now  |-D  functions  plotted  over 
all  possible  directions.  The  d  peaks  of  this  1-D  func¬ 
tions  constitute  a  good  estimate  of  the  d  unknown 
directions  (f,),4,,. 


I 
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Vin.  SIMULATION  RESULTS 

Because  of  space  limitations  we  will  present  only  the 
result  of  a  single  simulation.  The  array  consisted  of  Chree 
sensors  in  positions  (1.0).  (0.0)  and  (0.1)  arid  of‘  two 
sensors  in  positions  (-1.5)  and  (10.20).  The  signals 
were  independent  ARM  A  processes  with  identical  spectral 
density  centered  at  0.25  and  of  bandwidth  0.05.  The 
additive  noises  were  while  processes  and  the  signal-to- 
nolse  ratio  was  10  dB.  The  estimation  of  the  spectral 
density  matrix  was  done  by  the  periodogram  method 
using  an  FFT  of  64  points  with  a  total  of  300x04  samples. 
The  results,  using  the  estimator  given  by  (15a)  are  shown 
in  Fig.l.  The  two  peaks  corresponding  to  the  two  sources 
are  dearly  seen. 

DC  CONCLUDING  REMARKS 

A  new  non-parametric  approach  for  passive  localiza¬ 
tion  and  identiffcation  of  wideband  sources,  has  been 
presented.  The  approach  is  capable  of  handling  multiple 
correlated  sources  so  that  multipath  propagation  is 
included  as  a  special  case.  Dircction-of-ar rival  of  multi¬ 
ple  wideband  sources  is  also  included  as  a  special  case. 

The  method  presented  is  based  on  (he  cigenstruc- 
ture  of  the  spectral  density  matrix  of  the  received  sig¬ 
nals.  It  is  an  extension  of  the  method  of  Schmidt  (1979) 
and  Bienvenu  (1979)  to  the  case  of  wideband  sources.  It 
was  shown  that  the  eigenvector  sub-space  corresponding 
to  the  repeated  smallest  eigenvalue  of  the  spectral  den¬ 
sity  matrix  contains  all  the'  information  on  the  sources 
locations.  Algorithms  for  extracting  this  information 
from  an  estimate  of  the  spectral  density  matrix  of  the 
received  signals,  that  enable  trade-offs  between  resolu¬ 
tion  and  accuracy,  have  been  presented.  Simulation 
results  that  demonstrate  the  performance  of  the  pro¬ 
posed  algorithms  have  also  been  presented. 

We  should  note  that  though  the  proposed  method  is 
non-parametric.  in  the  sense  that  it  was  derived  for  non- 
parametrized  signals,  the  estimation  of  the  spectral  den¬ 
sity  matrix  of  the  received  signals  can  be  done  by  any 
method,  non-parametric  or  parametric,  whichever  is 
believed  to  give  a  better  estimate. 

Computation-wise,  the  method  is  quite  expensive 
since  it  involves  eigenvalue-eigenvector  decomposition  of 
the  spectral  density  matrix  at  each  frequency  bin  of  the 
received  signals.  However,  since  all  these  eigndecompo- 
sitions  can  be  carried  out  in  parallel,  using  special- 
purpose  hardware,  the  computation  time  can  be  reduced 
up  to  nearly  real-time. 
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